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Models of non-interacting fermions coupled to auxilliary classical degrees of freedom are relevant 
to the understanding of a wide variety of problems in many body physics, e.g. the description of 
manganites, diluted magnetic semiconductors or strongly interacting electrons on lattices. Monte 
Carlo sampling over the classical fields is a powerful, yet notoriously challenging, method for this 
class of problems - it requires the solution of the fermion problem for each classical field config¬ 
uration. Conventional Monte Carlo methods minimally utilize the information content of these 
solutions by extracting single temperature properties. We present a flat-histogram Monte Carlo al¬ 
gorithm that simulates a novel statistical ensemble which allows to acquire the full thermodynamic 
information, i.e. the partition function at all temperatures, of sampled classical configurations. 

PACS numbers: 05.10.Ln, 02.70.Ss 


Introduction - Bilinear Hamiltonians of lattice 
fermions coupled to classical degrees of freedom (contin¬ 
uous or discrete) are ubiquituous in contemporary many- 
body physics. These often arise as suitable approxima¬ 
tions in the description of systems where many different 
degrees of freedom contrive to yield complex and inter¬ 
esting physics. In these cases, some subsystem is treated 
classically as e.g. localized spins in double-exchange 
models [1K3J, models of Mn-doped (III, IV) semicon¬ 
ductors fj], the Ising t-J model @, adiabatic phonons 
in polaron models [fj-lsj, or one species of fermions in 
the Falicov-Kimball model fill ITTI| . Exactly solvable mod¬ 
els can also take this form, such as the seminal honey¬ 
comb lattice Kitaev model [l2| where interacting spins 
are mapped onto Majorana fermions coupled to static 
gauge fields. More generally, auxiliary field methods, e.g. 
based on the Hubbard-Stratonovich transformation, al¬ 
low to decouple fermion-fermion or fermion-boson inter¬ 
actions - the fields are then treated classically in con¬ 
junction with th e ap plication of the Suzuki-Trotter de¬ 
composition 13, fl4l | (see [l5| for a recent approximate 
scheme with static fields). 

The simplicity of the form of such models belies their 
complexity. Although obtaining eigenstates for fixed 
configurations of classical fields is computationally easy, 
summing over the exponential number of such configura¬ 
tions to obtain thermodynamic properties is notoriously 
expensive. An obvious approach towards this problem 
is via Monte Carlo (MC) sampling B]. Summing over 
the fermion states for fixed classical field configuration 
yields the conditional (grand) partition function. This 
quantity, at fixed temperature, serves as the Metropolis 
weight Jj| for performing a random walk in the space 
of classical field configurations []1|. The serious bottle¬ 
neck in these simulations is the repeated performance of 
the fermionic trace - and so exact diagonalization ED of 
the free fermion system - for executing the walk. Hence 


various improvements have been proposed to optimize 
the reevaluation of the weight - moment expansion of 
the fermion density of states by Chebyshev polynomials 
BEL low-rank matrix updates [21( or Green’s func¬ 
tions [221 ] Chebyshev expansion. On the other hand it 
seems naturally essential to optimize the extraction of 
information at each MC step. Indeed, while the expen¬ 
sive ED yields the conditional partition function at all 
temperatures, these are completely discarded by using 
only single temperature data in the above approaches. 
Here, we introduce an algorithm that fully exploits the 
thermodynamic information available at each diagonal¬ 
ization step in MC simulations of free fermions coupled 
to classical fields (FCCF). 


The paradigmatic Metropolis algorithm [17] suffers 
from critical slowing down at continuous phase transi¬ 
tions and prolonged trapping in metastable states at dis¬ 
continuous transitions. These can be overcome to a large 
extent using cluster update schemes [iij-Efi] or sampling 
extended ensembles such as Wang-Landau sam¬ 

pling of the density of states |32| 0- While these ap¬ 
proaches have been used in the study of classical and 
quantum systems, FCCF hold their own system-specific 
challenges rendering such applications difficult. In par¬ 
ticular, the effective Hamiltonian (corresponding to the 
energies in the Metropolis weight) of the classical fields 
in general contains temperature dependent, long-range, 
multi-particle interactions. Often, molecular dynamics 
(MD) or hybrid Monte Carlo methods using Langevin’s 
equation [Il-Eiil are better suited than standard MC 
simulations of such Hamiltonians although the accep¬ 
tance rate in these simulations crucially depends on the 
quality of the approximate action. Instead, we sample 
a novel extended ensemble bringing the advantages of 
Wang-Landau-like sampling to FCCF. 

The problem and method - We first generalize our con¬ 
siderations: a system is bi-separable if it may be sepa- 
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rated into two subsystems A and B such that for a given 
state of A all states of the system B can be efficiently 
summed over to obtain the conditional (grand) partition 
function Z(f3\A) or equivalently the conditional free en¬ 
ergy (grand potential) F(/3\A). This definition covers 
both FCCF (where subsystem A is classical) and many 
classical models (as e.g. the Ising model on a bipartite 
lattice where A and B are the spins on the two sublat¬ 
tices respectively). For bi-separable systems, the parti¬ 
tion function is decomposed as: 

Z(f3) = Trexp(-/3 H) = ^ (^ exp(-/3Ff B | A )) = 

A B 

^Z(/3|A)=^exp(-/3^(/3|A)) (1) 

A A 

i.e. the partition function Z(f3) is obtained by averag¬ 
ing the conditional partition functions over all configu¬ 
rations attainable in the system A. Notice that once 
the conditional energy spectrum is obtained, complete 
thermodynamic information associated with the expo¬ 
nentially large number of configurations of B encoded in 
exp(— j3F(/3\A)) for all inverse temperatures /3 becomes 
potentially available at each simulation step. 

Our problem can be stated as follows. Standard MC 
sampling for FCCF systems consists of walking, at fixed 
temperature, between different configurations of A with 
with Metropolis weights exp(—j3F(j3\A)). Here, we aim 
to obtain the entire partition function 0 by acquiring, 
at a given simulation step, the conditional partition func¬ 
tion Z(/3\A) for all arguments /? from the information 
(full spectrum) abundantly available for each fixed con¬ 
figuration on A. In principle this may also be achieved by 
parallel tempering (to efficiently sample configurational 
space of A) and a reweighting procedure, which however 
requires well chosen set of temperature intervals. Instead, 
in our method, we perform a random walk directly in the 
configurational space of A without referring to any spe¬ 
cific temperature. The basic challenge here is to obtain 
an appropriate importance sampling scheme over the (ex¬ 
ponential number of) configurations of A. 

The key behind any thermodynamic Monte Carlo sim¬ 
ulation lies in sweeping through energy space efficiently. 
A simple observation reveals that, for most systems, only 
a few configurations of A will lead to the spectrum of sys¬ 
tem energies containing the ground state energy. These 
configurations of A obviously must be effectively found 
by the importance sampling scheme. On the other hand 
the conditional DOSes pa{E ) of energy spectra associ¬ 
ated with typical configurations of A are expected to 
differ most also in their lower range. Therefore a key 
discriminant for configurations of A is the minimal en¬ 
ergy attainable by system at a given configuration which 
we will denote as £ m in(A). We consider two configura¬ 
tions of A as belonging to the same class when they have 
equal values of e m in- 


Formally, any importance sampling scheme can be rep¬ 
resented by assigning a weight function w(A) to the set 
of configurations of subsystem A, such that 


Z(f3) = '£w(A) 

A 


exp(-0F(0\A)) \ 

w(A) J ‘ 


( 2 ) 


The principle that we identify and implement for ac¬ 
quiring the partition function is that all minimal energy 
classes be visited by the algorithm. For this, notice that 
the minimal attainable energies e m i n can be associated to 
a density of states - p(e m in), which enumerates the num¬ 
ber of configurations of subsystem A attaining a given 
minimal energy. We emphasize that this ” auxiliary” DOS 
is distinct from the true DOS of the system, and does not 
determine the latter. 

Our algorithm consists of two separate sampling stages 
associated with first generating the weight distribution 
for importance sampling and then utilizing it to acquire 
thermodynamic information about the system, (i) The 
auxiliary DOS p(e m i n ) is readily obtained by performing 
a Wang-Landau (WL) algorithm [33] in the space of min¬ 
imal energies, where e m \ n plays the role of ’’energy” of a 
given configuration of A. (ii) Next, one performs a ran¬ 
dom walk in the space of configurations of A with weight 
function 


w(A) = l/p(e min (A)) (3) 

to sample Z{j.3) (Eq.©) for all (3. This choice of weight 
function not only allows to visit all classes of configura¬ 
tions but visits each class approximately the same num¬ 
ber of times yielding a flat histogram. 

The above principle can be viewed as a minimal nec¬ 
essary requirement for sampling energy space. It yields a 
coarse-grained view of energy space disregarding any dif¬ 
ferentiation between configurations belonging to a given 
class. Below, we show that the application of this prin¬ 
ciple leads to remarkably accurate results. 

Testing the algorithm - The convergence properties 
of the WL algorithm which constitutes the first stage of 
our algorithm have been widely studied. We note, that 
the random walk in the second stage of our algorithm 
(with fixed auxiliary DOS) fulfills, by fiat, detailed bal¬ 
ance. We test the convergence of our algorithm on the 
(a) Ising and (b) Potts model. Finally, we present contin¬ 
uous temperature results for the Falicov-Kimball model 
as a prominent example of FCCF. 

The Ising model with nearest neighbour interactions 
on a square lattice is an ideal benchmark for testing new 
algorithms since both an exact solution in the thermody¬ 
namic limit as well as large system high precision Monte 
Carlo results are available. This model is bi-separable 
which allows for direct application of our algorithm. In¬ 
deed for given configuration of spins on one of the sub¬ 
lattices (subsystem A) the conditional partition function 
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FIG. 1. (Color Online). The specific heat Cy per site for the 
2D Ising model obtained with our (black line) and WL (red 
line) algorithms. Lattice size is 60x60 with PBC. Inset shows 
the relative difference A Cy = (Cy i — Cy 2 )/Cyi. 

corresponding to all configurations of spins from the sec¬ 
ond sublattice (subsystem B) is simply: 

Z(P\A) = 2 Ar / 2 (cosh(2/3J)) JV2 (cosh(4/3J)) Ar4 , (4) 

where J is the coupling between spins, N /2 is the num¬ 
ber of sites on sublattice B. N 2 and N 4 are the number 
of spins on sublattice B subject to the nearest neighbour 
fields with absolute value 2J and 4J respectively. Of 
course, N 2 and IV4 are dictated by the configuration of 
spins on sublattice A. Fig. [j] shows the continuous tem¬ 
perature dependence of the specific heat, which is the 
temperature derivative of the free energy, obtained with 
our and Wang-Landau algorithms. These are generated 
for a 60 x 60 site lattice (with periodic boundary condi¬ 
tions PBC). The two curves show good agreement, indi¬ 
cating sufficiency of our prescribed importance sampling 
principle. Any essential problems regarding the effective¬ 
ness of this sampling scheme would have been expected 
to be apparent for this moderately large lattice size. 

We now highlight some important distinctions between 
the WL and our algorithms. A single step in the second 
stage of our algorithm entails accumulation of full ther¬ 
modynamic information from all configurations of sub¬ 
system B (a remarkable 2 1800 configurations in the pre¬ 
sented simulation) during each update move on subsys¬ 
tem A. In contrast the standard WL algorithm updates 
information from one configuration per move. However 
this ’’exponential update” comes at the initial cost of first 
obtaining the ’’auxiliary” DOS via a WL procedure. In¬ 
terestingly, the ’’auxiliary” DOS has to be determined 
with essentially a higher histogram flatness requirement 
than the system DOS in the direct WL simulation of the 
Ising model. This has two sources: (i) The ’’subsystem 
Hamiltonian” with energies e m in contains multi-spin and 
long range interactions unlike the original Ising interac¬ 
tion Hamiltonian, (ii) Sensitivity to the precision of the 


FIG. 2. (Color Online). The specific heat Cy per site for the 
2D 10-state Potts model obtained with our (black line) and 
WL (red line) algorithms. Lattice size is 30x30 with PBC. 
Upper inset - tips of the Cy curves from main panel. Lower 
inset - relative difference A Cy = (Cyi — Cy 2 )/Cyi. 

’’auxiliary” DOS. We have found that adding small fluc¬ 
tuations to /o(e m i n ) rapidly deteriorates results. Hence, 
our algorithm, should not be viewed as an alternative 
to the standard WL algorithm for classical models. Fi¬ 
nally while the WL algorithm directly outputs the system 
DOS, from which the partition function may be easily cal¬ 
culated, our algorithm yields the partition function. This 
difference has important practical consequences. In the 
WL algorithm the DOS values (which generally grow ex¬ 
ponentially with system size) are accumulated in a given 
run by multiplying them with the so-called modification 
factor / every time a given energy is visited. In subse¬ 
quent runs / is gradually reduced to 1. This ingenious 
trick allows to build up very big DOS values in a rea¬ 
sonable number of steps). In our procedure we sum the 
accumulated quantities exp(— /3F(/3\A))/w(A), and the 
values of both factors are already very big or very small. 
Therefore care must be taken to avoid roundoff error acu- 
mulation during summation. 

To illustrate that our algorithm succesfully simulates 
discontinuous phase transitions, we consider the 10-state 
Potts model on the square lattice with nearest neighbour 
interaction. Bi-separability here may be shown in similar 
fashion as in the Ising model. However the form of the 
conditional partition function is more complicated than 
© due to the multitude of values of local fields set by 
configurations of nearest neighbour Potts spins. Com¬ 
parison in Fig. [2] of the specific heat, obtained by the 
WL and our algorithm, for a 30 x 30 - site lattice with 
PBC reveals the effectiveness of our sampling scheme in 
this case as well. 

Now, consider the Falicov-Kimball Hamiltonian (FKH) 
of FCCF: 

HfK = -t 4' C j + U n i‘ n i ~ ^ n i- ’ 

<i,j> i i 
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FIG. 3. (Color Online). Specific heat Cy per site for the 2D Falicov-Kimball model at half-filling obtained with our (black line) 
and Metropolis (red circles) algorithms. The electron-ion onsite repulsion U/t = 0.25 (weak-coupling), 1 (moderate-coupling) 
and 8 (strong-coupling) from the left to the right panel. Lattice size is 16x16 with PBC. Insets: comparison of internal energies. 


where c\ and d\ are creation operators of mobile and 
immobile fermions respectively, = cjc,;, nf = djd,;. t is 
the nearest neighbour hopping integral, U is the Hubbard 
on-site interaction and fi is the chemical potential for c 
fermions (the number N& of immobile fermions is fixed). 
For a given configuration of immobile fermions {nf}, the 
set of single particle eigenenergies {ej} of mobile fermions 
is readily obtained, rendering the model bi-separable. 

Unlike the classical models above, an efficient standard 
WL algorithm is not directly applicable to the simula¬ 
tion of the FKH. We consider the FKH on a TV = L 2 
square lattice under PBC with L = 16, at half filling 
(fj, = U/2,Nd = N/ 2), for different values of U. Under 
these conditions, the FK model undergoes, at low tem¬ 
peratures, a transition to the charge density wave (CDW) 
ordered state, with Q = (7 r, 7 r) ordering wave-vector. 
The transition is continuous for large U with some evi¬ 
dence pointing to a change to discontinous transitions for 
small U 0, |4G ] . The application of our algorithm yields 
all-temperature results unlike standard Metropolis sam¬ 
pling over the immobile fermions. Moreover, since the 
Metropolis algorithm suffers from slow kinetics at dis¬ 
continous phase transitions our method should be par- 
ticulary useful in further investigations of the small U 
regime. Here, we restrict to proving that our sampling 
principle works. In Fig. [3l a comparison of results for the 
specific heat and internal energy obtained by Metropolis 
sampling and our algorithm is presented. On all diagrams 
the results are in good agreement within the accuracy of 
the local update Metropolis algorithm, which again in¬ 
dicates the effectiveness of our sampling. We mention 
here, that in obtaining the ” auxiliary” DOS in the first 
stage of our algorithm, we discretized values of minimal 
energies of the total system (but we have not discretized 
single particle energies anywhere else) by assigning them 
to (energy) bins of width 0.005< or O.OOlt and checked 
convergence. 

Summary and conclusions - In this Letter, we have 
presented a new Monte Carlo algorithm for problems of 


fermions coupled to classical degrees of freedom. This al¬ 
gorithm is based on Wang-Landau-like sampling in con¬ 
tradistinction to the commonly used Metropolis sam¬ 
pling. This allows in principle to overcome drawbacks 
of Metropolis schemes and importantly to fully exploit 
all available thermodynamic information at each diago- 
nalization step information that is mostly wasted in 
other MC schemes. The scheme is based on the no¬ 
tion of minimal energy attainable for a given classical 
field configuration. As a minimal requirement, we de¬ 
vised a rule that all such minimal energies be visited 
by the algorithm. This was achieved by first obtaining 
an auxiliary DOS for minimal energies via Wang-Landau 
sampling, and then using its inverse as the weight func¬ 
tion to perform a random walk in classical field config¬ 
urations accumulating the full partition functions. We 
benchmarked our recipe for several paradigmatic models 
of statistical physics achieving excellent agreement with 
known results. We mention here that while our principle 
of sampling minimal energies yields satisfactory results 
in the presented examples, more generally, supplemen¬ 
tary conditions may need to be identified in other cases. 
However, our results show that in principle accumulation 
of conditional partition functions at all temperatures at 
once using simple, temperature independent importance 
sampling is possible. An intriguing possibility is the ap¬ 
plication of this method to problems of quenched disor¬ 
der, where free energies need to be efficiently averaged 
over the disorder realizations. 


Finally, we comment on the most distinctive feature of 
our algorithm i. e. accumulation of full conditional par¬ 
tition functions per update. Standard algorithms update 
information from only single configurations per move. 
However, this is by no means the only possibility. N-fold 
way 41] algorithms may be seen as updating information 
from M configurations during each move (where M is the 
number of system sites). Our algorithm is remarkable in 
that it uniquely allows the update of information from an 
exponential number of configurations during each move 
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(~ exp(aM) for some constant a). Clearly, no algo¬ 
rithms exist that can update more information from a 
complexity point of view. 
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